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Particle size distributions in atmospheric clouds 

By Roberto Paoli & Karim Shariff f 


In this note, we derive a transport equation for a spatially integrated distribution function 
of particles size that is suitable for sparse particle systems, such as in atmospheric clouds. 
This is done by integrating a Boltzmann equation for a (local) distribution function over 
an arbitrary but finite volume. A methodology for evolving the moments of the integrated 
distribution is presented. These moments can be either tracked for a finite number of 
discrete populations (“clusters”) or treated as continuum variables. 


1. Introduction 

Particles are present in atmospheric clouds in several forms such as liquid droplets, non- 
volatile aerosols or ice crystals. Their microphysical properties control many processes 
such as the production of rain in stratocumulus clouds and radiation through cirrus 
clouds. These properties depend on the way particles interact with the surrounding air, 
through fluid-dynamic and thermodynamic processes. As these processes usually take 
place at small spatial scales, the interaction of particles with atmospheric turbulence is an 
important, though complex, problem in cloud physics (Shaw 2003). From a computational 
point of view, two major factors contribute to this complexity. First is the very high 
turbulence Reynolds number and the large range of spatial scales (Vaillancourt & Yau 
2000; Shaw 2003): for convective clouds, the ratio of energy-containing to dissipative 
length scales is 0(1O 5 ), while the Reynolds number of the largest eddies is 0(1O 6 to 
10'). The second factor is that the mean distance A between particles is of the order of 
the Kolmogorov scale 77 . Thus, if one contemplated direct numerical simulation (DNS) 
where all spatial scales are resolved, then one would have to track individual particles. 
Since DNS resolution is not affordable for these flows, Eulerian formulations for the 
liquid/solid phase are widely used in the simulation of clouds. These formulations fall 
into two main classes. The first is a “two fluid model” where particles are modeled as 
a continuum having a local mass density per unit volume. This approach carries no 
information about the distribution of the particle size. In the second approach, some 
physical properties (e.g. the mean radius) of some “ensemble” of particles are explicitly 
computed at each physical location x. The concept of particle size distribution at a 
point at this stage of our discussion is ambiguous but will be clarified later. A standard 
procedure used in two-phase flow models (e.g. Williams 1965; Cotton & Anthes 1989; 
Crowe et al. 1998) to describe an ensemble of particles is to define a distribution function 
/, in a manner analogous to the kinetic theory of gases. In kinetic theory, a distribution 
function /(x,v;f) is defined where /(x, v; t) <5x 5v represents the number of molecules 
that at time t. are between x and x + <5x and whose velocity is between v and v + <5v. 
It is assumed that, as (Sx, Sv ) — > 0, the phase volume still contains a sufficently large 
population of molecules that statistics can be used. This is usually true in gasdynamics 
because the mean free path of molecules is much smaller than the continuum scale one 
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cares about. The extension of this appraoch to particles other than molecules is formally 
straightforward (see for example the book by Williams 1965), as long as the continuum 
description remains valid. In the case of atmospheric clouds however, A ~ 77 so that only 
a few particles rather than a population are present in a volume V = 0(r] 3 ). The object 
of this note is first to derive a transport equation for an integrated distribution function 
To, describing an ensemble of particles inside an arbitrary but finite volume Vo . Then, 
this approach is specialized to atmopsheric clouds. Finally, a methodology is proposed 
to solve for the moments of the distribution function To- 


2. Distribution function in atmospheric clouds 

Let us consider a population of particles in a cloud from an ensemble of realizations. 
At any time t, each particle p occupies the position x p (t) in physical space, moves with 
velocity u p (t) and changes its radius r p (t). This population can be represented in a 
phase space, defined by the generalized coordinates q(f) and evolving via the generalized 
velocities U(q(i)): 


q (t) 


r{t) 

x(f) 

u {t) 


UMO) = § 


f(q(f)) 

u(q(f)) 

F(q (i)) 


(2.1) 


where x and u are the spatial coordinates and velocities, F is the functional law of the 
force per unit mass acting on the particle and r is the functional law of the growth rate 
of its radius. At any time t, each particle of the population occupies a point q(t) in this 
space. Instead of tracking each particle we wish to follow the evolution of a distribution 
function /(q(f); t) of the population. This is defined in such a way that 


/(q(f);t) 5Q{t) 


( 2 . 2 ) 


is the number of particles that at time t are inside a cube of volume SQ(t) in phase space, 
located between the coordinates q(f) and q(t) + hq(t). After a time dt, this volume has a 
value 5Q{t + dt) and the diagonally opposite corners of the cube are mapped to q (t + dt), 
and q(f + dt) + Sq(t + dt), respectively. At the same time, each particle can change its 
q and the particle number can vary because of evaporation or coagulation. To derive 
a transport equation for / we need to relate these quantities. For the sake of clarity, 
we will drop the explicit dependence on t in all variables, and define t' = t + dt, and 
q' = q(t + dt) . Then we have 

q' = q + U(q)df (2.3) 

q' + <Sq' = q+ (5q+ U(q + <5q) dt (2.4) 


The last term in (2.4) can be expanded as U(q + hq) = U(q) + (V q U)<5q where V q U is 
the gradient of U in phase space. Substituting (2.3) into (2.4) one gets 

<5q' = (I + dfV q U)(5q (2.5) 


The change of phase space volume SQ' — 6Q is related to the divergence of the generalized 
velocity U by 


SQ' - SQ 


SQ 


(V q • U) dt + 0(dt 2 ) 


( 2 . 6 ) 
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Let K (q; t) describe the general rate of gain or loss of particle number due to coagulation 
or evaporation, 


/(q';0 SQ' = /( q; t) SQ + A'(q; t)dt SQ 
Substituting (2.3), (2.4) and (2.6) into (2.7) and expanding gives: 


(2.7) 




V q /-U 


dt + 0(dt 2 ) 


1 + V q -U dt + 0(dt 2 ) 


SQ 


fSQ + KdtSQ (2.8) 


where all quantities are evaluated at t. Taking the limit dt — > 0 and neglecting higher 
order infinitesimals gives a Boltzmann equation 

^ + V q - (/U) = K (2.9) 

Finally, inserting the different components of q and U by means of (2.1) gives the more 
usual form employed in two-phase flow literature (e.g. Williams 1965): 

+ V x • (/u) + V u • (/F) + = K (2.10) 

The term F in (2.10) is the aerodynamic drag induced by the flow on particles, F ss 
(u p — u f)/r p (see Crowe et al. 1998) where u/ is the fluid velocity and t p = Ap p r 2 /\8p 
is a relaxation time. If the size of the particle r p is small, r p is also small and the particle 
velocity immediately adjusts to the flow velocity. In the following, we restrict our analysis 
to this case, so there is no dependence on F in (2.10): 

f+V x .(/u)+^ = A^ (2.11) 


3. Integrated distribution 

We now derive an integrated version of (2.11). Consider a point xo in physical space and 
an arbitrary (but finite) volume Vo( x o(t) around it. Then, a space integrated distribution 
function Tq{t, xo;f) can be defined as 


?o(r,xo;t) = [ f(r,x-,t)dV 0 (x) (3.1) 

Jv 0 

so that Ao(r, Xo; t) Sr represents the number of particles that at time t, are inside a finite 
volume Vo around xo and whose radius is in between r and r + Sr (the dimensions of this 
function are [SFq\ = A -1 whereas [/] = A -1 x A -3 ). The mean ijj 0 within the volume Vo 
of any quantity i p(r, x; t) associated with each particle is 


ip 0 (r,x 0 ;t) = 


f(r,x ; t)ip(r,x ; t) dV 0 (x) 


Tb(r,x 0 ;t) J Vo ’ 

Let us introduce a local coordinate y around xo, x = x(xo, y) = xo + y, so that 

V X • (•) = V XQ • (•) + Vy • (•). 

Using (3.2) for u and r, their integral values over Vo become 

no (r,x 0 ;t) = I f(r,x 0 +y;t)u(r,x 0 + y;t)dV 0 (y) 

Ao(r,x 0 ;t) J v 


/ f(r,x 0 +y-t)r(r,x 0 + y-t)dV 0 (y) 

JVo 


(3.2) 


(3.3) 

(3.4) 


r 0 (r, x 0 ; t) 


A 0 (r, x 0 ; t) 


(3.5) 
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Using (3.3) to express the divergence term V x • (/ u) = V Xo • (/ u) + V y • (/ u), a transport 
equation for To can be derived by integrating (2.11) inside volume Vjp 

J v d ~l t dV ° + y xo • [j v fudV^j + ^ Vy • (/ u) dV 0 + 

Ulj Mv °)=k Kdv ° <“> 

As in general xo and Vo vary in time, we need to switch volume integration and time 
derivative in the first term of the above equation. Using Leibnitz rule and (3.1), one has: 

(3j) 

where So is the surface around V 0 and S / is the velocity of S 0 with respect to the hxed 
reference frame. Using (3.7), Gauss theorem and the definitions (3.4) and (3.5) in (3.6), 
and introducing the relative velocity with respect to So, w = u — S/, finally gives 

+ V Xo • (T 0 u 0 ) + + j> f w • n 0 dS 0 = K 0 (3.8) 

which is a Boltzmann equation for the integral distribution function To- Equation (3.8) 
formally differs from (2.11) because of the surface integral in the left-hand side. This 
contains the (unknown) local distribution function /, which must be modeled in some 
way. Note, however, that if Vo is a material volume, then u = S / on So, and the surface 
flux goes to zero. 


Ergodic hypothesis 

In the previous derivation we had to introduce an ensemble of realizations in order to 
derive a local Boltzmann equation which we then integrated over a finite volume. We 
now make the hypothesis that Vo is large enough to contain a population of particles 
that we can by-pass the ensemble. In other words, we hypothesize that the integrated 
Boltzmann equation (3.8) is valid for a single realization if Vo is large enough. Typically, 
the grid size in cloud codes is A > 1 m while A ~ rj « 10 -3 m, so each grid cell contains 
at least 10 9 particles. 


3.1. Ensemble averages 

The total number of particles N 0 in spatial volume Vo can be obtained by integrating JF 0 
over all possible radii, 


/»oo 

A7 0 (x 0 ;t)= / , x 0 ; t) dr . 

Jo 


(3.9) 


Using (3.9), the ensemble average (ipo) of any variable ij) is obtained by integrating (3.2) 
over r : 

(^o) (x 0 ;t) = , / ,, / V’o(uxoU).7 r o(uxo;i)dr (3.10) 

A 0 (x 0 ; t) Jo 

In particular, the ensemble velocity and radius growth rate are 

1 r°° 

(uo) (x 0 ; t) = No ( Xo . t j J u o(u x 0 ; t) To (r, x 0 ; t) dr, 

(r 0 ) (x 0 ; t) = J ro(r, x 0 ; t) T 0 {r, x 0 ; t ) dr (3.11) 
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For further analysis, it will be useful to introduce the mean radius, (ro), and the variance, 
(Atq), of the population. Using (3.2), these are given by 

/»00 

(r 0 ) (x 0 ; t) = Nq ^ q .^ J rf 0 (r,xo\t)dr (3.12) 

/»oo 

( Arg) (x 0 ; t) = Nq (1- 0 . ^ J o l r - i r o) (u 0 ; t )] 2 T 0 (r, x 0 ; t) dr (3.13) 

Note that u(r, y; t) and (u 0 ) (x 0 ; t ) (same for r and (r 0 )) have a different physical mean- 
ing: u(y) represents the velocity of a particle in a neighboorood of y. On the other hand, 

Uo(xo) represents a statistical average within a population of particles and is a continuum 
velocity field, associated to any point, xo, of the physical domain. Indeed, one could, in 
principle, obtain (u 0 ) and (r 0 ) as 

, N 0 .N 0 

( u o) = W^' Up, ^ = N~ ^ ( 3 - 14 ) 

0 P = l 0 P = l 

where u p and f p are the velocity and the radius growth rate of particle p inside Vo- In 
the limit of N 0 —* oo, (3.11) and (3.14) are equivalent but we only have access to Tq 
because the details concerning u p and r p of each physical particle are unknown. 

The next step is to relate the continuum fields (uo) and (ro) to the corresponding flow 
variables. As we do not consider here any force acting on particles (Sec. 2), they are 
simply convected by the fluid. Therefore, there is no reason why two particles of the 
same population and different radius should have different velocities, i.e. Uo is statistically 
uncorrelated with r, 

u 0 (r, x 0 ; t) = (u 0 ) (x 0 ; t) = u/(x 0 ; t) (3.15) 

where u/(xo; t) is the fluid velocity at xo- The same arguments, cannot be applied to ro, 
i.e. fo yf (r) 0 , because each particle of the population may have a different growth rate 
due to different ’’reactions” to turbulent fluctuations in the flow-field, as discussed next. 


4. Particle growth by condensation 

The growth of the radius of a single particle in a medium at rest can be simply derived 
by considering a diffusion equation for water vapor on a particle surface (Pruppacher & 
Klett (1997) p. 502) and is given by 

dr = D (p v — p s v (T)) = PS 
dt r p w r r p w r 

where D is the diffusion coefficient of water vapor in air and T is the psychrometric 
correction associated with the latent heat of condensation; and p. w is the density of ei- 
ther water or ice. The vapor densities p v and p s v {T) are evaluated, respectively, at some 
“ambient” condition far from the particle and at the surface of the particle (which co- 
incides with the saturation value because vapor there is in thermodynamic equilibrium 
with water/ice). Thus, in such a single-particle picture, the radius growth rate is only 
controlled by the supersaturation S = p v — p s v . As first pointed out by Srivastava (1989) 
(see also Khvorostyanov & Curry 1999), this description cannot be extended straightfor- 
wardly to a population of particles. In fact, even in the absence of turbulence and uniform 
S initially, the available vapor in a cloud is not equally distributed among all particles 
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because of their random spatial distribution, so that the effective supersaturation avail- 
able at a droplet surface (“microsaturation”), can differ significantly from the overall 
ensemble averaged supersaturation (“macrosaturation”) (Srivastava 1989). In addition, 
in a turbulent cloud, each particle “reacts” in a different way to turbulent fluctuations in 
the flow- field: for example, if a supersaturation fluctuation arises, it will be absorbed by 
each particle through a complex diffusional process of vapor involving all elements of the 
population (Srivastava 1989). Several approaches have been developed in the atmospheric 
science literature (Pruppacher & Klett 1997) to try to solve this complex problem. One 
of these, the so-called “stochastic condensation” approach performs Reynolds averaging 
on the equation for condensational growth, resulting in covariances that can be though 
as “Reynolds stresses” (Shaw 2003). In particular, we follow Khvorostyanov & Curry 
(1999) (see also Pruppacher & Klett (1997) p.505) who use kinetic theory to relate the 
micro- and macro-saturation in a cloud. Their arguments are as follows. For the moment 
consider the situation where the particle radius is so small to be comparable with the 
mean free path of vapor molecules. In this case, one should account for the Brownian 
motion of molecules, that is the diffusion associated to the (random) molecular impact on 
particles surface. As shown by Pruppacher & Klett (1997), this can be done in (4.1) by 
introducing a modified diffusion coefficient D*(r ) which depends linearly on the radius r 
(Pruppacher & Klett (1997) p.506). It can be argued (Crowe et al. 1998; Khvorostyanov 
& Curry 1999) that the effects of turbulent fluctuations of vapor density or in super- 
saturation is similar to Brownian motion, whereas the molecular impact on particles is 
substituted by their interaction with turbulent eddies (note that this picture can also 
be extended to account for equivalent Brownian dispersion of particles, induced by fluc- 
tuating fluid forces rather than density fluctuations) (Crowe et al. 1998). The “micro” 
supersaturation S p available to particle p of an ensemble is 


S p = 


( So) 

(ro) 


(4.2) 


where (Sq) (x 0 ; t ) is the ensemble supersaturation available to the population within 
volume Vo- It represents the supersaturation that would be at xo if there were no particles, 
then it can be thught as the fluid supersaturation at xo, (So) (xo;t) = S/(xo;£). Using 
(4.2) and the previous formalism (r p — > r; S p — > S(r, y); f p — > r(r, y)), one gets to 


Substituting 

(ro) = 


S(r, y) = (So) = jS/y 
r (r 0 ) (r 0 ) 

f(r,y)= DS{r ' y) = DSf 
Tp w r Tp w (r 0 ) 


(4.3) and (4.4) into (3.11) finally gives 


7V 0 (x 0 ;£) Jo 


To ro dr 



DS(r, y) 
I 


dVo dr = 


DS f 
I >U, (ro) 


(4.3) 

(4.4) 


(4.5) 


5. Method of moments 

Even neglecting the surface term, (3.8) is a p.d.e. in four-dimensional space (r, xo;£) 
that can only be solved numerically. Some atmospheric cloud codes solve a transport 
equation for a distribution function by discretizing the particle size r in a finite number of 
bins, at each grid location (although it is not explicitly mentioned, they are conceptually 
discretizing (3.8) with neglected surface terms). 
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In this section, we present a simulation strategy based on the method of moments 
proposed in Paoli et al. (2002). The moments of the distribution To are defined by 


POO 

m fc (x 0 ;t)= / r k T 0 (r,x 0 -,t)dr 
Jo 


(5.1) 


Multiplying (3.8) by r k gives 

^ ( jf r k T o dr^j + V X0 • r k T 0 u 0 d?j = -T ( 0 r 0 r fe | + k J r k ~ l T 0 f 0 dr (5.2) 


Using (3.15), (4.5) and (5.1) and assuming that To — * 0 sufficiently fast as r — » oo, (5.2) 
becomes 

dnik _ . . DSt 

■^+V„-(u /mi ,)= — (5.3) 

Under all assumptions made, (5.3) describes the evolution of the moments of the integral 
distribution function To- An attractive property, deriving from the microsaturation model 
(4.2), is that the evolution of the k th moment only depends on the previous order moment 
which allows one to close the system (5.3) without any further assumptions and without 
presuming the shape of To- Using (3.9)-(3.13) and (5.1), the zero and the first two 
moments are easily found and are related to the ensemble average radius and variance, 


m 0 = N 0 , mi = N 0 (r 0 ) , 


m 2 


N 0 


(Arg) + (r 0 ) 2 . 


(5.4) 


The corresponding evolution equations are (we put a = D/Tp w ) 


+ V Xo • (u/ Nq) = 0 

dm\ , Nq 

— + V».(u ,m 1 ) = aS,- 

r\ 

+V Xo -(u / m 2 ) = 2a5 / lV 0 


(5.5) 

(5.6) 

(5.7) 


These equation are coupled to the continuum fluid phase through u f and Sf. In partic- 
ular, an increase in particle radius by condensation implies vapor depletion p v , with 

f°° A 2. -r , s f r , o N 0 m 2 /ko ^ 

Pv = - / 47T p w r r 0 To dr = — — / r T 0 dr = -47t p w abf (5.8) 

Jo T Pw (?’o) Jo ml 

The usual convection-diffusion-reaction equation for a scalar Y v in conservative form 
(p v = Y v pf where pf is the total gas phase density) then becomes: 


dp v 

~dt 


+ ^x 0 


(*U Pv) d^xo ' (p/TdVxpl^j) — p v 


—4:irp w aSf 


fV 0 m 2 

mi 


(5.9) 


Under all approximations made, (5.5)-(5.7) and (5.9) (together with Navier-Stokes equa- 
tions) describe the evolution of the first moments of the size distribution of a population 
of particles. These moments can be solved by using either an Eulerian or a Lagrangian 
description as discussed next. 


Eulerian description 

In this case, one has to solve for rrik using (5.5)-(5.7) with the further condition that the 
volume Vo is constant in time (in discretized form it can be the volume of a grid cell) . It is 
worth mentioning that To and mk are continuous functions of space (not grid averages!), 
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so their gradients contain all spatial fluctuations in a turbulent flow. In particular, if they 
are filtered in a LES approach, the correlations mpuj and m^Sf exist at subgrid scale 
level and must be modeled. 


An approximate Lagrangian method 

Let us devide the total number of particles in the cloud into N c ’’clusters” , each containing 
a fixed number, Nj, of particles (where j = 1 ,..,N C ). The position of each cluster is 
assumed to evolve according to 

^=“/( 4 ) ( 5 - 10 ) 

where Xg is the center of the volume Vg containing cluster j and u/(xq) is the fluid 
velocity at Xq. Note that we are assuming that the cluster advects rigidly without de- 
forming. Introducing the total derivative d()/dt = d()/dt + V Xo ()u/ in (5.5)-(5.7), one 
can “track” the moments of each cluster j as (note that the zeroth moment equation, 
Nj = const is now trivial) 


dm[ 

dt 


dm 3 2 

dt 



= 2 aS f Nj 


(5.11) 

(5.12) 


where S'/ = S/(x g) is the fluid supersaturation at Xg. The advantage of the Lagrangian 
approach is that the surface term in (3.8), which was neglected to get to (5.5)-(5.7), is 
now zero because Vg is a material volume, u = Sj- on Sg for all j. 


6. Conclusions 

In this note we derived a transport equation for the radius distribution function of a 
population of particles in an atmospheric cloud. We used a simple stochastic condensation 
model for the radius growth (taken from the atmospheric science literature) to relate the 
microsaturation around each particle to the macrosaturation of the entire population. 
Finally, we described a procedure to solve for the moments of the distribution, and showed 
that this can be either used in Eulerian continuum formulations or Lagrangian tracking 
of “clusters” of particles. Future plans include testing of this method by comparison 
with DNS of homogeneus and isotropic turbulence and individual particle tracking; and 
application of the method to natural and contrail-generated cirrus. 
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